package edu.gmu.atelier;


/**
 * Performs Eigen factorization (a.k.a. decomposition) on the passed matrix.
 * The method uses the power method described in Olver and
 * Shakiban Chapter 10.6.  Only finds the "dominant" eigen value/vector, i.e.
 * the one that dominates as exponent k -> inf.
 * 
 * @author James H. Pope
 */
public class EigenPowerFactor implements EigenFactor
{
    private Matrix a  = null;
    
    private Matrix d  = null;
    
    private Matrix s  = null;
    
    private Element eigenValue = null; // dominant eigen value
    private Vector  eigenVector= null; // vector for value, normalized
    
    public EigenPowerFactor( Matrix a )
    {
        this.a = a;
        
        /*
         * When to stop? When diff between eigen values is "small"
         */
        int maxIterations = 100;
        int i = 0;
        Epsilon e = Epsilon.makeDecimal(15);
        
        Vector vi = a.getFactory().makeVector( a.getRows() );
        vi.setElement(i, vi.getElement(0).getOne());
        eigenVector = vi;
        for (int k = 0; k < maxIterations; k++)
        {
            Vector nextEigenVector = a.mult(eigenVector);
            // Determine ratio of growth for (arbitrary) first component
            Element e0 = eigenVector.getElement(0);
            Element e1 = nextEigenVector.getElement(0);
            Element nextvalue= e1.mult(e0.inverse());
            //System.out.println("k="+k+" l="+eigenValue +" v="+eigenVector);
            // If we have "converged" bail, note test case 10.6.1 (e) bailed
            // after only one iteration, so this convergance method
            // not very good :(
            if( k > 1 )
            {
                double diff = eigenValue.sub(nextvalue).abs().doubleValue();
                if( e.similar(diff, 0.0) )
                {
                    System.out.println("Converged after "+k+" iterations.");
                    break;
                }
            }
            eigenValue = nextvalue;
            // normalize, value not changed, less chance of overflow
            eigenVector = nextEigenVector.normalize();
        }
    }
    
    public Matrix getA()
    {
        return this.a;
    }
    
    public Matrix getD()
    {
        return this.d;
    }
    
    public Matrix getS()
    {
        return this.s;
    }
    
    public Element getEigenValue()
    {
        return this.eigenValue;
    }
    
    public Vector getEigenVector()
    {
        return this.eigenVector;
    }
    
}
